subroutine usdefn( un, u,  ub, uf, ucb, ucsb, ucf, ucsf, pn,p, pb, pf, pcb, pcsb, pcf, pcsf, &
ron,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf,ds1,ds2,esm,dt,dtmin,epsilon2)
use Data_cfoam, only: epsilonPanikovski, epsMaximumPrinciple
implicit none
integer k,rank

real(8), dimension(3), intent(in) :: u, ub, uf, ucb, ucsb, ucf, ucsf
real(8), dimension(3), intent(out) :: un, esm

real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf
real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf
real *8 ds1,ds2,dt,dtmin
real *8 ut(3,9)
real *8 ei,z1,vtmp,dttmp
real *8 char3,q,r,qn,rn,u1x,ron,vf,wf,vb,wb
real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwellf,qmaxf,qminf
real *8 char3cf
real *8 char3csf
real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rmaxb,rminb
real *8 char3cb
real *8 char3csb
! SK EXTRA VARIABLES
real *8 qb,qcb,qcsb,rf,rcf,rcsf,rwellb
real *8 qwellb,qmaxb,qminb,rwellf,rmaxf,rminf
real *8 gl,gr,rnr,rnl,qnl,qnr,char1f,char1b,snd,sndf,sndb
real *8 smachl,smachr,smach,gn,rc1,rc2,rc3,rc4,qc1,qc2,qc3,qc4,var1,var2
! END SK
real *8 charc1,charc2,charc3
real *8 a(3,3)
real *8 epsilon
real *8 es(3)
real *8 deff,defb
real *8 dtmint,s,scb,scf,scsb,scsf,sn,sound
integer j
external deff
external defb
real *8 epsilon1,epsilon2
real(8) :: alfa, sn1,sn2, ut12,ut13,ut22,ut23, pn1,pn2, ron1,ron2

real(8), parameter :: gam = 1.4d0 ! heat capacity ratio for a diatomic gas
real(8), parameter :: gam2 = gam + gam, gamx = 2d0 * dsqrt( gam ) / ( gam - 1d0 ), gamy = ( gam - 1d0) / gam2   

! SK
! 1-u 2-ub 3-uf 4-ucb 5-ucsb 6-ucf 7-ucsf
call alcf_coor( esm, es, a )

ut( 1,1 ) = a( 1,1 ) * u(1) + a( 2,1 ) * u(2) + a( 3,1 ) * u(3)
ut( 2,1 ) = a( 1,2 ) * u(1) + a( 2,2 ) * u(2) + a( 3,2 ) * u(3)
ut( 3,1 ) = a( 1,3 ) * u(1) + a( 2,3 ) * u(2) + a( 3,3 ) * u(3)

ut( 1,2 ) = a( 1,1 ) * ub(1) + a( 2,1 ) * ub(2) + a( 3,1 ) * ub(3)
ut( 2,2 ) = a( 1,2 ) * ub(1) + a( 2,2 ) * ub(2) + a( 3,2 ) * ub(3)
ut( 3,2 ) = a( 1,3 ) * ub(1) + a( 2,3 ) * ub(2) + a( 3,3 ) * ub(3)

ut( 1,3 ) = a( 1,1 ) * uf(1) + a( 2,1 ) * uf(2) + a( 3,1 ) * uf(3)
ut( 2,3 ) = a( 1,2 ) * uf(1) + a( 2,2 ) * uf(2) + a( 3,2 ) * uf(3)
ut( 3,3 ) = a( 1,3 ) * uf(1) + a( 2,3 ) * uf(2) + a( 3,3 ) * uf(3)

ut( 1,4 ) = a( 1,1 ) * ucb(1) + a( 2,1 ) * ucb(2) + a( 3,1 ) * ucb(3)
ut( 2,4 ) = a( 1,2 ) * ucb(1) + a( 2,2 ) * ucb(2) + a( 3,2 ) * ucb(3)
ut( 3,4 ) = a( 1,3 ) * ucb(1) + a( 2,3 ) * ucb(2) + a( 3,3 ) * ucb(3)

ut( 1,5 ) = a( 1,1 ) * ucsb(1) + a( 2,1 ) * ucsb(2) + a( 3,1 ) * ucsb(3)
ut( 2,5 ) = a( 1,2 ) * ucsb(1) + a( 2,2 ) * ucsb(2) + a( 3,2 ) * ucsb(3)
ut( 3,5 ) = a( 1,3 ) * ucsb(1) + a( 2,3 ) * ucsb(2) + a( 3,3 ) * ucsb(3)

ut( 1,6 ) = a( 1,1 ) * ucf(1) + a( 2,1 ) * ucf(2) + a( 3,1 ) * ucf(3)
ut( 2,6 ) = a( 1,2 ) * ucf(1) + a( 2,2 ) * ucf(2) + a( 3,2 ) * ucf(3)
ut( 3,6 ) = a( 1,3 ) * ucf(1) + a( 2,3 ) * ucf(2) + a( 3,3 ) * ucf(3)

ut( 1,7 ) = a( 1,1 ) * ucsf(1) + a( 2,1 ) * ucsf(2) + a( 3,1 ) * ucsf(3)
ut( 2,7 ) = a( 1,2 ) * ucsf(1) + a( 2,2 ) * ucsf(2) + a( 3,2 ) * ucsf(3)
ut( 3,7 ) = a( 1,3 ) * ucsf(1) + a( 2,3 ) * ucsf(2) + a( 3,3 ) * ucsf(3)

char3   = ut(1,1)
char3b  = ut(1,2)
char3f  = ut(1,3)
char3cb = ut(1,4) 
char3csb= ut(1,5)
char3cf = ut(1,6) 
char3csf= ut(1,7)

s    = dlog(p)    - gam * dlog(rho)
sb   = dlog(pb)   - gam * dlog(rhob)
sf   = dlog(pf)   - gam * dlog(rhof)
scb  = dlog(pcb)  - gam * dlog(rhocb)
scsb = dlog(pcsb) - gam * dlog(rhocsb)
scf  = dlog(pcf)  - gam * dlog(rhocf)
scsf = dlog(pcsf) - gam * dlog(rhocsf)

gf = gamx * dexp( scsf / gam2 )
gb = gamx * dexp( scsb / gam2)

soundf = dsqrt( gam * pcsf / rhocsf )
soundb = dsqrt( gam * pcsb / rhocsb )
sndb   = dsqrt( gam * pb   / rhob )
snd    = dsqrt( gam * p    / rho )
sndf   = dsqrt( gam * pf   / rhof )


! SK USE PANIKOVSKI FOR THE CENTRED EXPANSION FAN 
! RIGHT CELL
epsilon = epsilon2
if(char3.le.0 .and. char3f.gt.0)epsilon = epsilonPanikovski !0.7 ! empirical calculation methods for important points ! SK
if(char3+snd.le.0 .and. char3f+sndf.gt.0)epsilon = epsilonPanikovski !=0.7
if(char3-snd.le.0 .and. char3f-sndf.gt.0)epsilon = epsilonPanikovski !=0.7

! LEFT CELL
epsilon1 = epsilon2
if(char3.ge.0 .and. char3b.lt.0)epsilon1 = epsilonPanikovski !=0.7 ! if(char3b.le.0 .and. char3.gt.0)epsilon1=0.7
if(char3b+sndb.le.0 .and. char3+snd.gt.0)epsilon1 = epsilonPanikovski !=0.7
if(char3b-sndb.le.0 .and. char3-snd.gt.0)epsilon1 = epsilonPanikovski !=0.7
	
! END CENTRED FAN -- THIS SHOULD BE PLACED BEFORE ANY SUBSEQUENT REDEFINITION OF CHAR3F, ETC

q    = char3    - gf * p   **gamy
qf   = char3f   - gf * pf  **gamy
qcf  = char3cf  - gf * pcf **gamy
qcsf = char3csf - gf * pcsf**gamy

char2f=char3csf-soundf	
qwellf=(qcsf-qcf)/(0.5d0*dt)+char2f*(qf-q)/ds2
qmaxf=dmax1(q,qcf,qf)+dt*qwellf + epsMaximumPrinciple
qminf=dmin1(q,qcf,qf)+dt*qwellf - epsMaximumPrinciple

qnl=(2d0*qcsf-qf*(1.0d0-epsilon))/(1.0d0+epsilon)
if(qnl>qmaxf) qnl=qmaxf
if(qnl<qminf) qnl=qminf

!SK  + 4 TERMS IN THE OTHER DIRECTION  (B)
q    = char3    - gb * p   **gamy      
qb   = char3b   - gb * pb  **gamy
qcb  = char3cb  - gb * pcb **gamy
qcsb = char3csb - gb * pcsb**gamy

! SK + EXTRAPOLATION AND CORRECTION IN THE OTHER DIRECTION (B)
char2b=char3csb-soundb	
qwellb=(qcsb-qcb)/(0.5d0*dt)+char2b*(q-qb)/ds1
qmaxb=dmax1(q,qcb,qb)+dt*qwellb + epsMaximumPrinciple
qminb=dmin1(q,qcb,qb)+dt*qwellb - epsMaximumPrinciple
qnr=(2d0*qcsb-qb*(1.0d0-epsilon1))/(1.0d0+epsilon1)
if(qnr>qmaxb) qnr=qmaxb
if(qnr<qminb) qnr=qminb

! SK
r    = char3    + gb * p   **gamy
rb   = char3b   + gb * pb  **gamy
rcb  = char3cb  + gb * pcb **gamy
rcsb = char3csb + gb * pcsb**gamy

char1b=char3csb+soundb
rwellb=(rcsb-rcb)/(0.5d0*dt)+char1b*(r-rb)/ds1	
rmaxb=dmax1(r,rcb,rb)+dt*rwellb + epsMaximumPrinciple
rminb=dmin1(r,rcb,rb)+dt*rwellb - epsMaximumPrinciple

rnr=(2d0*rcsb-rb*(1.0d0-epsilon1))/(1.0d0+epsilon1)
if(rnr>rmaxb) rnr=rmaxb
if(rnr<rminb) rnr=rminb

!SK +4 TERMS IN THE OTHER DIRECTION (F)
r    = char3    + gf * p   **gamy
rf   = char3f   + gf * pf  **gamy
rcf  = char3cf  + gf * pcf **gamy
rcsf = char3csf + gf * pcsf**gamy

! SK -- THE SAME FOR THE OTHER DIRECTION (F)

char1f=char3csf+soundf
rwellf=(rcsf-rcf)/(0.5d0*dt)+char1f*(rf-r)/ds2	
rmaxf=dmax1(r,rcf,rf)+dt*rwellf + epsMaximumPrinciple
rminf=dmin1(r,rcf,rf)+dt*rwellf - epsMaximumPrinciple
rnl=(2d0*rcsf-rf*(1.0d0-epsilon))/(1.0d0+epsilon)
if(rnl>rmaxf) rnl=rmaxf
if(rnl<rminf) rnl=rminf

    ! SK ALL CHAR'S AT THE RHS BELOW MUST BE REFERRED TO THE CENTRES 
    sound=soundf+soundb
    charc3=char3csf+char3csb
    charc1=char3+sound   
    charc2=char3-sound



    ! SK DEFINITION //////////////////////////////////////
    gr=gb
    gl=gf
    char3f=char3csf
    char3b=char3csb
    ! END DEFINITION//////////////////////////////////////

    ! SK THE BASE CASE: SUBSONIC LEFT/RIGHT GOING WAVES//////////////////////

    gf=gl
    gb=gr
    rn=rnr
    qn=qnl


! SK MINIMUM VARIATION -- CHAR1


    if(char1f<0.0 .and. char1b>=0.0) then
        rc1=rcb-dt*char1b*rnl/ds1
        rc2=rcf+dt*char1f*rnl/ds2
        var1=dabs(rc1-rc2)

        rc3=rcb-dt*char1b*rnr/ds1
        rc4=rcf+dt*char1f*rnr/ds2
        var2=dabs(rc3-rc4)

        if(var1.lt.var2) then
            rn=rnl
            gf=gl
        else
            rn=rnr
            gf=gr
        end if
    end if

 ! SK END OF MINIMUM VAR CHAR1 

! SK MINIMUM VARIATION -- CHAR2


    if(char2f<0.0 .and. char2b>=0.0) then
        !	QCSLN1=QCL-DT*QLN*CHAR2L/SMESHL*DLNNL
        !	QCSN1=QC+DT*QLN*CHAR2R/SMESH*DLNN
        qc1=qcb-dt*char2b*qnl/ds1
        qc2=qcf+dt*char2f*qnl/ds2
        var1=dabs(qc1-qc2)

        !   QCSLN2=QCL-DT*QRN*CHAR2L/SMESHL*DLNNL
        !	QCSN2=QC+DT*QRN*CHAR2R/SMESH*DLNN
        qc3=qcb-dt*char2b*qnr/ds1
        qc4=qcf+dt*char2f*qnr/ds2
        var2=dabs(qc3-qc4)

        if(var1.lt.var2) then
            qn=qnl
            gb=gl
        else
            qn=qnr
            gb=gr
        end if
    end if

 ! SK END OF MINIMUM VAR CHAR2 

    u1x=qn+gf*(rn-qn)/(gf+gb)
    call GetAlfa(u1x,alfa)
    if(alfa>=1.0d0) then   ! if(charc3<=0.0) then ! SK changed, 2011-05-17
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
        ut(1,8)=0.0d0
        ut(2,8)=deff(ut(2,1),ut(2,3),ut(2,6),ut(2,7),char3csf,dt,ds2,epsilon)
        ut(3,8)=deff(ut(3,1),ut(3,3),ut(3,6),ut(3,7),char3csf,dt,ds2,epsilon)
        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
        ron=(pn/exp(sn))**(1.0d0/gam)
    else if(alfa<=0.0d0) then
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        ut(1,8)=0.0d0
        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon1)
        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon1)
        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
    else
        sn1 = deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
        ut12 = deff(ut(2,1),ut(2,3),ut(2,6),ut(2,7),char3csf,dt,ds2,epsilon)
        ut13 = deff(ut(3,1),ut(3,3),ut(3,6),ut(3,7),char3csf,dt,ds2,epsilon)
        pn1 = ((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
        ron1 = (pn1/exp(sn1))**(1.0d0/gam)

        sn2 = defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        ut22 = defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon1)
        ut23 = defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon1)
        pn2 = ((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
        ron2 = (pn2/exp(sn2))**(1d0/gam)

        sn = alfa * sn1 + (1.0d0 - alfa) * sn2
        ut(1,8)=0.0d0
        ut(2,8) = alfa * ut12 + (1.0d0 - alfa) * ut22
        ut(3,8) = alfa * ut13 + (1.0d0 - alfa) * ut23
        pn = alfa * pn1 + (1.0d0 - alfa) * pn2
        ron = alfa * ron1 + (1.0d0 - alfa) * ron2
    end if
!    if(u1x<=0.0) then   ! if(charc3<=0.0) then ! SK changed, 2011-05-17
!        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
!        ut(1,8)=0.0d0
!        ut(2,8)=deff(ut(2,1),ut(2,3),ut(2,6),ut(2,7),char3csf,dt,ds2,epsilon)
!        ut(3,8)=deff(ut(3,1),ut(3,3),ut(3,6),ut(3,7),char3csf,dt,ds2,epsilon)
!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
!        ron=(pn/exp(sn))**(1.0d0/gam)
!        !if (k.eq.5051) write(14,*) 'f',pn,ron
!        !if (k.eq.5051) write(14,*) 'fqr',q,r
!    else 
!        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
!        ut(1,8)=0.0d0
!        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon1)
!        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon1)
!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
!        ron=(pn/exp(sn))**(1d0/gam)
!        !if (k.eq.5051) write(14,*) 'b',pn,ron
!        !if (k.eq.5051) write(14,*) 'bqr',q,r
!    end if

! SK END OF SUBSONIC ////////////////////////////////////////////

! SK RIGHT-GOING SUPERSONIC/////////////////////////////////////////////


    if(char1b>0.0 .and. char2f>0.0) then
        gb=gr
        gf=gb
        rn=rnr
        qn=qnr

        u1x=0.5d0*(rn+qn)
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        ut(1,8)=0.0
        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon1)
        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon1)
        pn=((rn-qn)/2.0d0/gb)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
    end if


! SK END OF R/G SUPERSONIC ////////////////////////////////

! SK LEFT-GOING SUPERSONIC/////////////////////////////////////////////

    if(char1b<=0.0 .and. char2f<0.0) then
        gb=gl
        gf=gb
        rn=rnl
        qn=qnl

        u1x=0.5d0*(rn+qn)
                sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
	        ut(1,8)=0.0
	        ut(2,8)=deff(ut(2,1),ut(2,3),ut(2,6),ut(2,7),char3csf,dt,ds2,epsilon)
	        ut(3,8)=deff(ut(3,1),ut(3,3),ut(3,6),ut(3,7),char3csf,dt,ds2,epsilon)
        pn=((rn-qn)/2.0d0/gf)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1.0d0/gam)
    end if

! SK END OF L/G SUPERSONIC ////////////////////////////////

! SK SPECIAL TREATMENT FOR THE SONIC POINT -- RIGHT-GOING WAVE

    if(char1b>0.0 .and. char2b<=0.0 .and. char2f>=0.0) then

!SMACHL=WICL/SOUNDL
!SMACH=WIC/SOUND
!SM=0.5*(SMACHL+SMACH)
!G=DSQRT(GAM*DEXP(SRN)**(1./GAM))
!PNN=(RRN/(AL+SM*G))**(2.*GAM/(GAM-1.))
!UNN=SM*G*RRN/(AL+SM*G)
!RONN=(PNN/DEXP(SRN))**(1./GAM)

        gb=gr
        gf=gb
        rn=rnr
        qn=qnr
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        ut(1,8)=0.0
        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon1)
        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon1)
        gn=dsqrt(gam*dexp(sn)**(1./gam))

        smachl=char3b/soundb
        smachr=char3f/soundf
        smach=0.5*(smachl+smachr)

        pn=(rn/(gb+gn*smach))**(2.0d0*gam/(gam-1.0))
        u1x=smach*gn*rn/(gb+gn*smach)
        ron=(pn/exp(sn))**(1.0d0/gam)


    end if

! SK END OF SONIC POINT -- R/G WAVE ////////////////////////////////

! SK SPECIAL TREATMENT FOR THE SONIC POINT -- LEFT-GOING WAVE

    if(char1b<=0.0 .and. char1f>=0.0 .and. char2f<0.0) then

        gb=gl
        gf=gb
        rn=rnl
        qn=qnl
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
        ut(1,8)=0.0d0
	        ut(2,8)=deff(ut(2,1),ut(2,3),ut(2,6),ut(2,7),char3csf,dt,ds2,epsilon)
	        ut(3,8)=deff(ut(3,1),ut(3,3),ut(3,6),ut(3,7),char3csf,dt,ds2,epsilon)
        gn=dsqrt(gam*dexp(sn)**(1./gam))

        smachl=char3b/soundb
        smachr=char3f/soundf
        smach=0.5*(smachl+smachr)

        !G=DSQRT(GAM*DEXP(SLN)**(1./GAM))
        !PNN=(-QLN/(AR-SM*G))**(2.*GAM/(GAM-1.))
        !UNN=-SM*G*QLN/(AR-SM*G)

        pn=(-qn/(gb-gn*smach))**(2.0d0*gam/(gam-1.0))
        u1x=-smach*gn*qn/(gb-gn*smach)
        ron=(pn/exp(sn))**(1.0d0/gam)

    end if

! SK END OF SONIC POINT -- L/G WAVE ////////////////////////////////


call transformr(ut(1:3,8),ut(1:3,9),a)

!do j=1,3    
!    un(j)=ut(j,9)+u1x*es(j)
!end do

un(1) = ut(1,9) + u1x * es(1)
un(2) = ut(2,9) + u1x * es(2)
un(3) = ut(3,9) + u1x * es(3)

vtmp  = dsqrt( un(1) * un(1) + un(2) * un(2) + un(3) * un(3) )
sound = dsqrt( gam * pn / ron )

dtmint= min( ds1, ds2 ) / ( sound + vtmp )
if (dtmint.lt.dtmin) dtmin = dtmint

return 
end subroutine usdefn
